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The harmonic balance method (HBM) was originally developed for finding periodic solutions 
of electronical and mechanical systems under a periodic force, but has later been adapted 
to self-sustained musical instruments. Unlike time-domain methods, this frequency-domain 
method does not capture transients and so is not adapted for sound synthesis. However, 
its independence of time makes it very useful for studying every periodic solution of the 
model, whether stable or unstable without care of initial conditions. A computer program 
for solving general problems involving nonlinearly coupled exciter and resonator, "Harmbal" , 
has been developed based on the HBM. The method as well as convergence improvements 
and continuations facilities are thorougly presented and discussed in the present paper. Ap- 
plication of the method is demonstrated on various problems related to a common model 
of the clarinet: a reed modelled as a simple spring with and without mass and damping, a 
nonlinear coupling and a cubic simplification of it, and a cylindrical bore with or without 
dissipation and dispersion as well as a bore formed as a stepped cone. 

PACS numbers: 43.75.Pq, 43.58.Ta 



I. INTRODUCTION 

Since Helmholtz, 7 it has become natural to describe a 
self-sustainedi musical instrument as an exciter coupled 
to a resonator. More recently, Mclntyre et al. ? have 
highlighted that simple models are able to describe the 
main functioning of most self-sustained musical instru- 
ments. These models rely on few equations whose im- 
plementation is not CPU-demanding, mainly because the 
nonlinearity is spatially localized in an area small com- 
pared to the wavelength. This makes them well adapted 
for real-time computation (including both transient and 
steady states). These models are particularly popular in 
the framework of sound synthesis. 

On the other hand, calculation in the frequency do- 
main is suitable for determining periodic solutions of the 
model (the values of the harmonics as well as the playing 
frequency) for a given set of parameters. Such infor- 
mation can be provided by an iterative method named 
the harmonic balance method (HBM). Though the name 
"harmonic balance" seems to date back to 1936, 7 the 
method was popularized nearly forty years ago for electri- 
cal and mechanical engineering purposes, first for forced 
vibrations, 7 later for auto-oscillating systems. 7 The 
modern version was presented rather shortly after by 
Nakhla and Vlach. 7 In 1978, Schumacher was the first 
one to use the HBM for musical acoustics purposes with 
a focus on the clarinet. 7 However in this paper, the play- 
ing frequency is not determined by the HBM. This short- 



coming is the major improvement brought by Gilbert et 
al. 7 eleven years later, who proposed a full study of the 
clarinet including the playing frequency as an unknown 
of the problem. 

The fact that the HBM can only calculate periodic so- 
lutions, may seem as a drawback. Certainly, transients 
such as the attack are impossible to calculate, and the pe- 
riodic result is boring to listen to and does not represent 
the musicality of the instrument. Therefore the HBM 
is definitely not intended for sound synthesis. Neverthe- 
less, self-sustained musical instruments are usually used 
to generate harmonic sounds, which are periodic by def- 
inition. The HBM is thus very useful to investigate the 
behavior of a physical model of an instrument, depend- 
ing on its parameter values. This is possible for both 
stable and unstable solutions, without care of precise ini- 
tial conditions. Moreover, HBM results can be compared 
to approximate analytical calculations (like the variable 
truncation method (VTM)), 7 in order to check the va- 
lidity of the approximate model considered. 

The present paper is based on the work of Gilbert et 
al. 7 Our main contributions are: extension of the diver- 
sity of equations managed, improved convergence of the 
method, introduction of basic continuation facilities, and 
from a practical point of view, faster calculations. 

While the main idea is already described by Gilbert 
et al., 7 Section [H] details the principle of the HBM, in 
particular the discretization of the problem, both in time 
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and frequency. 

Section IIIII is devoted to the various contributions of 
the current work, which are applied in a computer pro- 
gram called Harmbal. 7 The framework is defined to in- 
clude models with three equations: two linear differential 
equations, written in the frequency domain, and a nonlin- 
ear coupling equation in the time domain (see Sec. If 11 A*j l. 
As usual in the HBM, this system of three equations is 
solved iteratively. The solving method chosen (Newton- 
Raphson, Sec. If 11 B"|) has been investigated and its conver- 
gence has been improved through a backtracking scheme 
rSecs. lTTTT1 a,nd lTTTT)t . 

To illustrate the advantages of the HBM and the im- 
provements, a few case studies were performed and are 
presented in Section IIVI They are based on a classi- 
cal model of single reed instruments which is presented 
in Section llVAl In Sections llV Bl and further, simplifica- 
tions to each of the three equations are introduced so that 
the results could be compared to analytical calculations, 
both for cylindrical and stepped-cones bores. Finally the 
full model is compared to time-domain simulations. This 
also shows the modularity of Harmbal. The comparison 
is achieved through the investigation of bifurcation dia- 
grams as the dimensionless blowing pressure is altered. 
The derivation of a branch of solution is obtained thanks 
to basic continuation with an auto- adapt at ive parameter 
step. 

Finally, various questions are tackled through practi- 
cal experience from using Harmbal. Section Ivl discusses 
multiplicity of solutions and poor robustness in the fre- 
quency estimation. 

II. NUMERICAL METHOD 

A. The harmonic balance method 

The harmonic balance method is a numerical method 
to calculate the steady-state spectrum of periodic solu- 
tions of a nonlinear dynamical system. In this paper we 
are only concerned with periodic solutions. The following 
provides a detailed and general description of the method 
for a nonlinearly coupled exciter-resonator system. 

Let X(uik), k = 0, . . . , N t — 1 be the Discrete Fourier 
transform (DFT) of one period x(t), < t < T, of a 
T-periodic solution of a mathematical system to be de- 
fined. X(o;fe) will have a number of complex components 
Nt, which depends on the sampling frequency f s = 1/T S 
with which we discretize x(t) into N t = T/T s equidistant 
samples. Furthermore, uJk—2Tif p T s k is the angular fre- 
quency of each harmonic of the fundamental frequency 
f p of the oscillation, referred to as the playing frequency. 
Note that the sampling frequency f s = N t f p is automati- 
cally adjusted to the current playing frequency so that we 
always consider one period of the oscillation while keep- 
ing Nt constant. Note also that Nt should be sufficiently 
large to avoid aliasing. Moreover, if it is chosen a power 
of two, the Fast Fourier transform (FFT) may be used. 
Assuming that N p < N t /2 harmonics is sufficient to de- 



scribe the solution, we define X € M 2A p+ 2 as the N p + 1 
first real components (denoted by 5ft) of X(oJk) followed 
by their imaginary components (5): 

A = [5R(A(^ )),...,5ft(A(^ p )), (I) 

Note that the components Xq and Ajv p +i are the real and 
imaginary DC components respectively (and that Xn +i 
is always zero). Our mathematical system can thus be 
defined by the nonlinear function F : R 2Ar p+ 3 — > R 2Ar p+ 2 : 

X = F(XJ P ). (2) 

Until now, the playing frequency has silently been as- 
sumed to be a known quantity. In autonomous systems, 
however, the frequency is an additional unknown, so that 
the iVp-harmonic solution seeked is defined by 2N p +3 un- 
knowns linked through the 2N p +2 equations How- 
ever, it is well known that as A is a periodic solution of 
a dynamical system, any X' deduced from A by a phase 
rotation (i.e. a shift in the time domain) is also a solution. 
Thus an additional constraint has to be added in order 
to select a single periodic solution among the infinity of 
phase-rotated solutions. A common choice (see Ref.? ) 
is to consider the solution for which the first harmonic 
is real (i.e. its imaginary part, A^r +2 , is zero). This ad- 
ditional constraint decreases the number of unknowns to 
2N p +2 for an A p -harmonic periodic solution. Thus we 
get F : R 2N + 2 -> R 2N + 2 , and it is now possible to find 
periodic solutions, if they exist. 

Finally, a simple way of avoiding trivial solutions to 
equation J3J) is to look for roots of the function G : 

R 2N p +2 ^ R 2JV p+2; defined by 

a(x-,u) = *- f ^' fp \ (3) 

i.e. G(X, f p ) = 0. This equation is usually solved numer- 
ically through an iteration process, for instance by the 
Newton- Raphson method as in our case. How to handle 
the playing frequency f p will be discussed in the following 
section. 

B. Iteration by Newton-Raphson 

The equation G{X,f) — 0, G being defined by equa- 
tion J3J), is nonlinear and has usually no analytical so- 
lution. (For readability we leave out the index p on the 
playing frequency until end of Sec. IIIII ) This section de- 
scribes the common, iterative Newton-Raphson method. 
This is the method used in the program Harmbal (see 
Section lHl|) although it had to be refined with a back- 
tracking procedure to improve its convergence, as dis- 
cussed in Section lill Dl 

For the sake of later reference, it is useful to recollect 
the principles of Newton's method for a one-dimensional 
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FIG. 1. The iteration process of Newton's method 



problem g(x) — 0. Starting with an estimate x° of the 
solution, the next estimate a; 1 is defined as the intersec- 
tion point between the tangent to g at xq and the a;-axis. 
The method can be summarized as 



x l+1 =x l ~ 



g'{x 1 )' 



(4) 



This is repeated, as shown in Figure ^ while increasing 
the iteration index i until g(x l ) < e, where e is a user- 
defined threshold value. 

In our 2A^,+2-dimensional case, we have a vector prob- 
lem: we search (X, f) for which G(X, f) — 0. Newton's 
method is generalized to the Newton-Raphson method, 
which may be written 7 : 

(x i+i j i+i ) = (x^n-iJhy'-dix^n (5) 

where 3 l G =\/G{X l , p) is the Jacobian matrix of G at 
{X % , P). Note that all derivatives by Xm + 2, which was 
chosen to be zero, are ignored. The column N p +2 in 
the Jacobian is thus replaced by the derivatives with re- 
spect to the playing frequency /. J Z G is thus a (2iV p +2)- 
square matrix. This means that line number N p +2 
in equation J5J gives the new frequency / instead of 
We define the Newton step AX=X l+1 -X l 



X 



N p +2- 



(where A/=/ l+ — /' replaces AXn p +2), which follows 
the local steepest descent direction. 

The Jacobian may be found analytically if G is given 
analytically, but it is usually sufficient to use the first- 
order approximation 



dGi 



jk 



G 3 {X + 5X k J)-G 3 {XJ) 



dX k SX 
except for k = N p + 2, in which case we use 



J 



j,N+2 



dGi „ GjjXJ + Sfi-GjjXJ) 
df Sf 



(6) 



(7) 



The components of 8X k are zero except for the fcth one, 
which is the tiny perturbation 5X. The iteration has 
converged when \G l \ = \G{X\ f)\ < e. We found e = 
1CP 5 to be a good compromize between computation time 
and solution accuracy. 



III. IMPLEMENTATION AND HARMBAL 

A. Equations for self-sustained musical 
instruments 

Though, to the authors' knowledge, the harmonic bal- 
ance method in the context of musical acoustics with un- 
known playing frequency has only been applied to study 
models of clarinet-like instruments, it should be possible 
to consider many different classes of self-sustained instru- 
ments. It is well accepted that sound production by a 
musical instrument results from the interaction between 
an exciter and a resonator through a nonlinear coupling. 
Moreover, in most playing conditions, linear modelling of 
both the exciter and the resonator is a good approxima- 
tion. 

Therefore, within these hypotheses, any musical in- 
strument could be modelled by the following three equa- 
tions: 



Z e (ui)X e (uj) = X c {uj) (a) 
X c (u) = Z r (uj)X r (uj) (b) 
F(x c (t),x e (t),x r (t)) = Q (c) 



(8) 



where Z e is the dynamic stiffness and Z r is the input 
impedance of the exciter and the resonator, respectively, 
and X e and X r are the spectra describing the dynamics of 
the exciter and the resonator during the steady state (pe- 
riodicity assumption) . X c is the spectrum of the coupling 
variable. All these quantities, and thus equations JHJvi 
b), are defined in the Fourier domain. Equation JSJi) is 
written in the time domain, where J- is & nonlinear func- 
tional of x c , x e , and x r , which are the inverse Fourier 
transforms of X c , X e , and X r , respectively. We apply 
the discretization as described in Section III Al implying 
that equations b) become vector equations where the 
impedances must be written as real (2iV p +2) x (2iV p +2)- 
matrices to accommodate the rules of complex multipli- 
cation: 



U> \Q(Z(f)) R(i(/)) J 



(9) 



where 



Z(f) 



( Z(p) 
Z(u)i) 

\ 



(10) 



Z{un p ) J 



is complex, and 3t(Z) and $s{Z) are the real and imag- 
inary components of Z. The system © is solved itera- 
tively by Harmbal according to the scheme illustrated in 
Figure 

In Harmbal, these equations are easily defined by writ- 
ing new C functions. Only superficial knowledge of the 
C language is necessary to do this. 

Three cases related to models of single reed instru- 
ments with cylindrical or stepped-conical bores are stud- 
ied in particular in Section IIVI in order to validate the 
code and to illustrate the modularity of Harmbal. 
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X c (t) 



Eq.©) 



DFT" 1 
X e ,f x e (t) 



x r (t) 
| DFT 

Xr,f 

| Eq.JHt) 



I AX C ,A/- 



N-R 



FIG. 2. The iteration loop of the harmonic balance method 
for a musical instrument (notations defined in the text) 



B. Practical characteristics of Harmbal 

Both fast calculation, good portability, and indepen- 
dence of commercial software are easily achieved by pro- 
gramming in C, whose compiler is freely available for 
most computer platforms. It is, however, somewhat dif- 
ficult to combine portability with easy usage, because an 
intuitive usage normally means a graphical and interac- 
tive user interface, while the handling of graphics varies 
a lot between the different platforms. 

We have chosen to write Harmbal with a nongraphical 
and non-interactive^ user interface. The major advantage 
of this is that independent user interfaces may be further 
developed depending on need. 

Our concept is to save both the parameters and the 
solution in a single file. This file also serves as input 
to Harmbal while individual parameters can be changed 
through start-up arguments. The solution provided by 
the file works as the initial condition for the harmonic 
balance method. Thus the lack of a simple user interface 
is compensated by a simple way of re-using an existing 
solution to solve the system for a slightly different set of 
parameters. Solutions for a range of a parameter values 
may thereby be calculated by changing the parameter 
stepwise and providing the previous solution as an ini- 
tial condition for the next run. The Perl script hbmap 
provides such zeroth-order continuation facilities. This 
procedure may also be used when searching for a solu- 
tion where it is difficult to provide a sufficiently good 
initial condition, for instance by successively increasing 
N p when wanting many harmonics. 

C. Convergence of Newton-Raphson 

When merely employing the Newton-Raphson method 
to determine the solution of the system at a given set 
of parameters, we have found that it is impossible to 
find a solution at particular combinations of the param- 
eters. Indeed, for the clarinet model of Section IIVB II 
no convergence was obtained for particular values of the 
parameter 7 (the dimcnsionlcss blowing pressure) and its 
neighborhood. This is seen as discontinuities, or holes, 
in the curves in Figure (see Section IIVI for the under- 
lying equations and parameters). Note that the solu- 




FIG. 3. Solution holes: first pressure harmonic Pi versus 
blowing pressure 7 for different N p with Nt = 128, £ = 0.5, 
and 77 = 10 -3 . (Even N p give the same as N p — 1.) Equations 
and parameters are defined in section ITvl 



Y= 0.4194 
Y= 0.4195 
T= 0.4196 
y= 0.4197 
y= 0.4198 



0.2096 0.2097 0.2098 0.2099 0.21 0.2101 0.2102 0.2103 0.2104 
Pi 

FIG. 4. Gi as Pi varies around the solution Gi — for 
various 7 around a hole at 7 ~ 0.4196. Nt = 128 and N p = 1. 



tions seem to go continuously through this hole and that 
the positions of the holes and their extent vary with the 
number of harmonics N p taken into account. The curves 
were calculated by the program hbmap. In this case we 
have decreased 7 from 0.5 downward in steps of 10 -4 
and drawn a line between them except across 7 values 
where solution failed. In the holes, the Newton-Raphson 
method did not converge, either by alternating between 
two values of P (i.e. X c ) or by starting to diverge. 

To study the problem, we simplified the system to a 
one-dimensional problem by setting N p = 1 , thus leaving 
Pi as the only nonzero value. G\ thus became the only 
contributor to |G|, and a simple graph of G\ around the 
solution G\ — could illustrate the problem, as shown 
in Figure 0] We see that the curve of G\(Pi) has in- 
flection points (visible as "soft steps" on the curve) at 
rather regular distances. At the centre of a convergence 
hole, i.e. for 7 ~ 0.4196, an inflection point is located at 
the intersection with the horizontal axis. This is a school 
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example of a situation where Newton's method does not 
converge because the Newton step APi brings us alter- 
natingly from one side of the solution to the other, but 
not closer. 

In fact, the existence of inflection points is linked with 
the digital sampling of the continuous signal. If the sam- 
pling rate is increased, i.e. if N t is increased, the steps 
become smaller but occur more frequently, as shown for 
N t = 32, 128, and 1024 in Figures EJi-c. The derivative 
dGi I dP\ is included in the figures to quantify the impor- 
tance of the steps. According to the Figurcs[SJi-c it seems 
reasonable to increase Nt to avoid convergence problems. 
However, this would significantly increase the computa- 
tional cost. Another solution is therefore suggested in 
the following. 

D. Backtracking 

When the Newton-Raphson scheme fails to converge, 
it often happens because the Newton step AX leads 
to a point where \G(X,f)\ is larger than in the pre- 
vious step. However, acknowledging that the Newton 
step points in the direction of the steepest descent, there 
must be a point along AX where \G(X,f)\ is smaller 
than in the previous iteration of the HBM. A backtrack- 
ing algorithm described in Numerical Recipes 7 (Sec. 9. 7) 
solves the problem elegantly by shortening the Newton 
step as described here. The principle is illustrated in the 
simple one-dimensional case in Figure El where g(x) re- 
places \G(X,f)\, although we use the multidimensional 
notation in the following. Defining the A axis along the 
Newton step, we simply take a step XAX in the same 
direction, where < A < 1. The optimal value for A is 
the one that minimizes the function h(X): 

h(X) = i|G(A l + AAA)| 2 (11) 

with derivative 

h 'W = { J o-G)\ jti+XAji .AX. (12) 

During the calculation of the failing Newton step, we 
computed G(X l ) and G(X l+1 ), so now it is possible 
to calculate with nearly no additional computational ef- 
fort h{0) = WGiX 1 )] 2 , h'(0) = -\G{X l )\ 2 , and h(l) = 
^G^+AX)] 2 = \\G{X i+1 )\ 2 . This allows to propose 
a quadratic approximation of h for A between and 1, 
for which the minimum is located at 



h(l)-h(0)-h'(0)' v ' 

It can be shown that Ai should not exceed 0.5, and in 
practice Ai > 0.1 is required to avoid a too short step at 
this stage. 

If \G{X l + XiAX)\ still is larger than \G(X% h(X) is 
then modelled as a cubic function (using h(X±) which 



has just been calculated). The minimum of this cu- 
bic function gives a new value A 2 , again restricted to 
O.lAi < A 2 < 0.5Ai. This calculation requires solving 
a system of two equations, so if also X 2 is not accepted 
because \G(X l + A2AX)| is still too large, we do not en- 
hance to a fourth-order model of h, which would increase 
the computational cost much more. Instead, subsequent 
cubic modellings are performed using the most two recent 
values of A. In practice, however, not many repetitions 
should be necessary before finding a better solution, if 
possible. 

IV. CASE STUDIES 

A. The equations for the clarinet 

The three equations (JSJi-c) may be constructed by 
physical modelling. In the case of the clarinet, a com- 
mon simple model is described below. We limit the de- 
scription in the following to a brief presentation based on 
dimensionless quantities, dimensional variables being de- 
noted by a hat (~) hereafter (see Fritz et al. ? for further 
details). 

The exciter is an oscillating reed which may be mod- 
elled as a spring with mass and damping: 

y + g e y + vly = —(p-Pm.), (14) 

Me 

where y is the dynamic reed displacement, and p and p m 
are the dynamic pressure in the mouthpiece, i.e. the in- 
ternal pressure, and the static blowing pressure in the 
player's mouth, respectively. The constants g e , and 
w e represent the mass per area, the damping factor, and 
the angular resonance frequency of the exciter (the reed). 
The dots over y denote the time derivative. In dimension- 
less form, equation 114(1 becomes 

Mx + Rx + Kx = p, (15) 

where p = p/pm and x = y/H + ^/K with 7 = p m /PM- 
The equilibrium reed opening is H as shown in Figure 
In the static regime, when blowing harder than a maxi- 
mum pressure pm, i-e. p m > Pm (7 > 1), the reed blocks 
the opening, i.e. y — —H, so we get p — and can con- 
clude that K — 1 for the current reed model. 

Like Fritz et al. ? we relate the dimensionless time to 
the resonance angular frequency uj r of the resonator (the 
bore), i.e. t = u) r t, so that the values of the dimensionless 
mass M, damping i?, and spring constant K become 

K = [i e Huj 2 /pM = l, (16) 
R = Kg e u r /u 2 = g e uj r /u 2 e , (17) 
M = KuPjul = (18) 

In the Fourier domain, Equation i|15|) thus takes the form 
of equation Z e (u>)X{u>) = P(oj), where 

W = l-MQ\iRQ, (19) 
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FIG. 6. The principle of backtracking in one dimension. Ob- 
jective is to estimate the root of g(x) (solid curve). Broken 
lines with arrows show how the Newton step Ax from x leads 
to divergence. h(X) (dot-dashed curve) is a 2nd order expan- 
sion of g(x) along the Newton step, i.e. the A axis. Minimum 
of h(X) should be closer to the root of g(x) than g(x + Ax). 




FIG. 7. Illustration of the mouthpiece 



for i = \f—l and ui = 2-klo /uj r — uj / ' f r is the dimension- 
less angular frequency in the Fourier domain. 

A common minimum model for the clarinet assumes a 
simple reed with no mass or damping, thus M = R = 0. 
Equation H15|) reduces to x = p. 

The resonator (i.e. the air column in the bore of the 
instrument) is commonly described by its frequency re- 
sponse Z r (uj). For a simple cylindrical bore of length I 
with a closed and an ideal open end, the resonance fre- 
quencies are odd multiples of f r — c/Al, c being the sound 
speed in the air column. 7 The input impedance of the 
bore may thus be expressed in dimensionless quantities 



as 

Z r (uj) = = i tan(^ + (1 - i)a(wj) , (20) 

where a(oj) = ipr]y [ /uj/2ir with ip ~ 1.3 for common con- 
ditions in air and r/ being the dimensionless loss param- 
eter, which depends on the tube length, typically 0.02 
for a normal clarinet with all holes closed. Zq = pej S is 
the characteristic impedance of the cylidrical resonator, 
S being its cross section, and p the density of air. The 
last term in the argument of equation (|20|l includes the 
dispersion as the real part and viscous losses as the imag- 
inary part. 

Equation JSJs) becomes 

P(lu) = Z r {uj)U(Lu), (21) 

where P(co) and U(lu) = U(lu)Zq/pm are the dimension- 
less internal pressure and volume flow of air through the 
mouthpiece in the Fourier domain. 

The coupling equation l|Ht), is given by the Bernoulli 
theorem with some supplementary hypotheses applied 
between the mouth and the outlet of the reed chan- 
nel. The coupling equation is nonlinear and must be 
calculated in the time domain. This leads to the follow- 
ing expression for the dimensionless airflow through the 
mouthpiece 7 : 

u(p, x) = C (1 + x - 7) yyy - P\ sign(7 - p) (22) 

as long as 1 > 7 — 1, and u = otherwise. £ = 
ZqwH ppM is a dimensionless embouchure parame- 
ter roughly describing the mouthpiece and the position 
of the player's mouth, w being the width of the open- 
ing and p the density of the air. C, is also related to the 
maximum volume velocity entering the tube. 7 

If the reed dynamics were not taken into account, we 
had x = p and thus 

u(p) = C(l+P-7)v / l7-p|sign(7-p) (23) 
for p > 7 — 1, and, as before, u — otherwise. 
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B. Verification of method and models 

In the following we want to verify that the HBM (and 
its implementation in Harmbal) gives correct results. By 
using very low losses in the resonator (small rj) we can 
compare the results of the HBM with analytical results. 
Rising the attenuation in the resonator and including 
mass and damping for the exciter, we compare with nu- 
merical results from real-time synthesis of the same sys- 
tem. This also gives us the opportunity to illustrate the 
modularity of Harmbal as we change the models of the 
resonator and the nonlinear coupling. 

1. Helmholtz oscillation for cylindrical tubes 

To compare the HBM results with analytical results, 
we assume a nondissipative, nondispersive air column, 
i.e. setting rj = and thus a = in equation (|20|l . 
Furthermore, we assume that the reed has neither mass 
nor damping and thus use equation (|23|l . The resulting 
square- wave amplitude (the Helmholtz motion) 7 may be 
found by solving u(p) = u(—p), which results from the 
fact that the internal pressure p(t) and the power p(t)u(t) 
averaged over a period are zero according to the lossless 
hypothesis. 7 This leads to the square oscillation with 
amplitude 



p(j) = V~37 2 + 4 7 - 1. 



(24) 



This result is compared with the results calculated by 
Harmbal (for the same set of equations, but r\ = 1CP 5 
instead of 77 = to avoid infinite impedance peaks) for 3, 
9, 49, and 299 harmonics close to the oscillation threshold 
in Figure [SJ and at 7 = 0.4 in Figure EO which is far from 
the threshold. 

As expected, the solution using the HBM shows good 
convergence towards the Helmholtz motion as the num- 
ber of harmonics increases. Note the deviation for higher 
harmonics close to the threshold, even for 299 harmon- 
ics. Dissipation in the resonator (77 = 10~ 5 7^ 0) causes 
higher harmonics to be damped more in this area of 7 
than for higher blowing pressures (as explained e.g. in 
Ref.? ). The deviation from a square- wave signal is thus 
more noticeable close to the threshold, and as the HBM 
calculations imposed a nonzero dissipation, this is prob- 
ably the reason for the small deviation in Figure |SJ The 
deviation is not visible in the time domain. 

A popular simplification of the nonlinear function 1 2. 'ill 
is a cubic expansion for small oscillations (e.g. Ref.) 7 7 7 : 



u(p) = u 00 + Ap + Bp 2 + Cp 3 , 



(25) 



where uqq, A, B, and C are easily found by expanding 
equation (|23(l . Its Helmholtz solution is easily calculated 
like above, yielding 



'8 7 2 (37- 1) 
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FIG. 8. The Helmholtz solution, eq. I12H compared with the 
HBM truncated to 3, 9, 49, and 299 harmonics close to the 
oscillation threshold (7 = 0.334, £ = 0.5, 77 = 10~ 5 ). (a) 
frequency domain, (b) time domain. 



The influence of the difference between the two versions 
of the nonlinear function is investigated in Figure 1101 for 
the lossless case. Close to the oscillation threshold, Fig- 
ure HOfa . we see that there is no significant difference 
between the two versions of the nonlinear equation, as 
expected. The fact that the HBM is lower for higher 
harmonics is as before due to the small attenuation we 
had to include to perform the numerical calculations. Far 
from the threshold, however, Figure HTIb. we see that the 
cubic expansion fails to approximate the nonlinear equa- 
tion. For lower harmonics this error is larger than the 
attenuation effect in the HBM calculations. This is fur- 
ther discussed by Fritz et al. 7 

In Figure ^] we have completed some of the curves 
that we failed to make in Figure|21 and even increased the 
number of harmonics, owing to the backtracking mecha- 
nism. Admittedly, at N p — 49, a few holes can still be 
seen, but the convergence is significantly improved. 

Here the amplitude of the first harmonic is plotted 
for different numbers of harmonics as a function of the 
blowing pressure 7 together with first harmonic of the 
Helmholtz solution, deduced from equation Q24[l. In prac- 
tice, the solution at 7 = 0.4 was found and then hbmap 
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FIG. 9. The Helmholtz solution, eq. 1241 compared with the 
HBM truncated to 3, 9, 49, and 299 harmonics far from the 
oscillation threshold (7 = 0.40, ( = 0.5, rj = 10" 5 ) in the 
frequency domain. 
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FIG. 10. The (lossless) Helmholtz motion and the (almost 
lossless) HBM for 299 harmonics using the full nonlinear- 
ity (12311 and the cubic expansion (1251 (a) close to the os- 
cillation threshold (7 = 0.334) and (b) far from it (7 = 0.40) 
for £ = 0.5, 77 = 10~ 5 . Above the 11th harmonic only every 
10th harmonic is shown. 



was used to make Harmbal calculate solution for each of a 
large number of subsequent values of 7 down to the oscil- 
lation threshold by using the previous solution as initial 
value. The procedure was repeated from 7 = 0.4 up to 
the point where the reed started to beat, i.e. for p < 7 — 1 
in equation l|23|l . Without losses (Helmholtz solution) 
the beating threshold does not arrive before 7 = 0.5, 
and this should be expected for the nearly lossless case 
studied with the HBM also. However, the number of har- 
monics N p taken into account in the HBM calculations is 
too small to follow the sharp edges of the square signal. 
The resulting overshoots in p(t), as seen in Figure |SJd, 
cause p to prematurely exceed the criterion for beating. 
The beating threshold converges to 0.5 as N p increases 
(see also Ref.? ). Note that, for the chosen value of 
it can be calculated following Hirschberg 7 (eq.(45)) that 
above 7 ~ 0.45, the Helmholtz solution loses its stabil- 
ity through a subharmonic bifurcation (a period-doubling 
occuring). 

By Figure UTI we can also verify that the model expe- 
riences a direct Hopf bifurcation (which is known since 
the work of Grand et al.). ? Thus, a single harmonic is 
enough to study the solution around the threshold. Far 
from the threshold, more harmonics have to be taken 
into account for Pi to converge toward the Helmholtz so- 
lution. This is not obvious and for example contradictory 
with the hypothesis made for the VTM. ? Thus Harmbal 
appears as an interesting tool to evaluate the relevance of 
approximate methods according to the parameter values. 

2. Helmholtz oscillation for a stepped conical tube 

The saxophone works similarly to the clarinet, but the 
bore has a conic form. In this section we compare the 
HBM calculations with analytical results, and in order to 
calculate the Helmholtz motion when losses are ignored, 
we need to simplify the cone by assuming that it consists 
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of a sequence of N sylinders of length I and cross section 
Si = + 1)51, Si = S being the cross section of the 
smallest cylinder, and i — 1,...,N (see Ref.? ). The 
total length of the instrument is thus L = Nl. The input 
impedance of such a stepped cone may be written as 

Zr M = cot(f -ia( W '))+cot(^-ia(£)) ' (2?) 

where a/ = 2w/(iV+l) when cj = 2irf/f r , where f r is the 
first eigcnfrcqucncy of this resonator. We have ignored 
the dispersion term here. Equation (|27[1 is used instead 
of equation Q2U|). and the damping a(cd) = i/jt) ^uj/2tt 
is zero in the analytic Hclmholtz case and very small 
(77 = 2 • 10~ 5 below which convergence became difficult) 
for the calculations with the HBM. 

As before, the pressure amplitude of the ideal lossless 
case is calculated by solving u(p) — u(—Np), and two 
solutions are possible: 3 



p±( 7 )- 



± 



(i y_l)(2-3 7 ) 

2( N 2 -N+l) 

y/(iV-l) 2 + (^+1) 2 (~3 7 2 +47^T) 

2(7V 2 -7V + 1) 



(28) 



as long as 7 < l/(N + 1) for the standard Helmholtz 
motion (p + ) and 7 < N/(N + 1) for the inverted one 
(p~), which is unstable. Above these limits p + = 7 and 
p~ = —j/N. The magnitude of the first harmonic of a 
square or rectangular wave is then given by 



sin 



N+l 



P ± {l)- 



(29) 



For N = 1, equation (|28|l reduces to equation (|24|l . For 
higher N, the pressure oscillation becomes asymmetric. 
We take the case N = 2 and get 



^(7) 



2 - 3 7 ± V-277 2 + 36 7 - 



(30) 



This result is compared with HBM calculations in Fig- 
ure El for 7 = 0.31. Theoretically, the spectrum of the 
Helmholtz solution, Figure 112b . shows that every third 
component is missing (actually zero) while the remain- 
ing components decrease in magnitude thus forming the 
asymmetric pressure oscillation as shown in Figure IT^k . 
The HBM, on the other hand, suggests that the first 
component in each pair be smaller than the second com- 
ponent. This results in a dip at the middle of the long, 
positive part of the period (i.e. on both extremities t = 
and t — 1024 of the curve in Figure EJ. The same was 
observed for N — 3 and N = 4, where the long part of 
the period was divided by similar dips into three and 
four parts, respectively (not shown). The number of 
time samples, N t did not change this fact, but as Fig- 
ure E| indicates, the dips gradually become narrower as 
the number of harmonics N p increases. This indicates 
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FIG. 13. Amplitude of first harmonic Pi as a function of 
the blowing pressure 7 for the Helmholtz solution 13UH for 
2-stepped cone and the HBM truncated to 2, 5,. . . , and 63 
harmonics, the last coinciding with Helmholtz (£ = 0.2, r\ = 
2 • 10 -5 ). Only nonbeating regimes are shown. 



J. Acoust. Soc. Am. 



S. Farner et al.: Harmonic balance calculations for musical instruments 9 



that the HBM approaches the Helmholtz solution as N p 
approaches infinity. 

A bifurcation diagram is plotted in Figure 1131 Sim- 
ilarly to Figure 1111 for the cylindrical bore, the ampli- 
tude of the first harmonic is plotted for different num- 
ber of harmonics as a function of the blowing pressure 
7. The Helmholtz solution (equation 1(23)1 with N = 2) 
is also plotted. As shown by Ollivier et al., ? the lower 
part of the upper branch and the branch of the inverted 
Helmholtz motion are unstable. 

In practice, these curves are more difficult to obtain 
with hbmap than for the cylindrical bore, especially close 
to the subcritical oscillation threshold around 7 = 0.28, 
where computation was not possible at this low losses. 
More sophisticated continuation schemes should be con- 
sidered to obtain complete curves. However, it is obvi- 
ous from the diagram that the model experiences a sub- 
critical Hopf bifurcation, which agrees with the conclu- 
sion of Grand et al.." This means that a single-harmonic 
approximation is not enough to study the solution around 
this threshold, since the small-amplitude hypothesis does 
not hold. Further from the threshold, convergence to- 
ward the Helmholtz motion is ensured as the number of 
harmonics N p is increased. 

Only the nonbeating reed regime is considered in the 
figure and, similarly to Figure ^2 it can be noted that 
the beating threshold for the model with N p harmon- 
ics depends on N p but converges toward the Helmholtz 
threshold 7 = 1/3 (corresponding to the lossless, contin- 
uous system) as N p is increased. 

3. Validation with time-domain model 

When adding a mass and damping to the reed or vis- 
cous losses and dispersion to the pipe, it is more diffi- 
cult to compare Harmbal results with analytic solutions. 
This has been done by Fritz et al. ? as far as the play- 
ing frequency is concerned, by comparison with approxi- 
mate analytical formula. Here, we propose to confront 
both the playing frequency and the amplitude of the 
first partial with numerical results obtained with a time- 
domain method. We use a newly developed (real-time) 
time-domain method (here called TDM) by Guillemain et 
al.. 7 It is based on the same set of equations as presented 
in Section II V Al except that the impedance of the bore is 
slightly modified to be expressed as an infinite impulse 
response. In the Fourier domain, it can be expressed as 



Z r (Cj) 



1 — die 



b e~ 



uD 



1 — a\e 



b e~ 



■jD ■ 



(31) 



where uj — Coj f s , f s being the sampling frequency, and 
the integer D = round(/ s /2/ r ) the time delay in samples 
for the sound wave to propagate to the end of the bore 
and back. The constants ai and bo are to be adjusted 
so that the two first peaks of resonance have the same 
amplitude as the two first peaks of equation (|20|l . 

To express equation (|31|) using our terminology, we 



TABLE I. The values of M and R for three strengths of 
reed interaction. The bore parameters are D — 247 (f r = 
103.4Hz), ai — 0.899, and bo — 0.0946 for sampling frequency 
f s = 51100 Hz. 



Reed 


cj e /Hz 




M 


R 


Weak 


10000 


0.1 


1.070-10" 4 


1.034- 10" 3 


Normal 


2500 


0.2 


1.712-10~ 3 


8.28- 10~ 3 



remember that lu = 2nf/f r and obtain 



Z t (lS) — 



1 — aie 



% e 



-iuj/2 



1 — a±e 



fr 

"^77 



5 e 



-iuj/2 



(32) 



In this section, we also include the mass and damp- 
ing of the reed, so M and R are no longer zero. The 
TDM does not work for M = R = 0, or even for values 
close to this, so we have used a reed with weak interac- 
tion with the pipe resonance as well as one with close to 
normal reed impedance. The corresponding values for Lo e 
and q e = g e /Lo e are shown in Table [I] Figure IT^k shows 
the bifurcation diagram for two values of £ and for weak 
and normal reed impedance, while Figure fHb shows the 
corresponding variation in the dimensionless playing fre- 
quency f p /f r - The lines represent the continuous solu- 
tions of the HBM, and the symbols show a set of results 
derived from the steady-state part of the TDM signal. 
The TDM symbols fall well on the lines of the HBM, ex- 
cept for £ = 0.50 when 7 approaches 0.5. Then the TDM 
experiences period doubling, i.e. two subsequent periods 
of the signal differ. At the same time, not being able to 
show subharmonics, the HBM shows signs of a beating 
reed, possibly a solution that is unstable and thus not 
attainable by time-domain methods. 

Note that three points have to be verified before com- 
paring results from the HBM and the TDM: 

The numerical scheme used in the TDM to approxi- 
mate the time derivatives in the reed equation Ijl5|l re- 
quires discretization. Depending on the sampling fre- 
quency f s , the peak of resonance of the reed deviates 
more or less from the one given by the continuous equa- 
tion. For normal reed interaction (/ e =2500 Hz), the de- 
viation is negligible, but it may become significant in 
the case of weak reed interaction, where the peak is at 
10000 Hz. However, the fact that the reed and the bore 
interact weakly in the latter case, implies that the exact 
position of the peak has little importance. Therefore, 
at the used sampling frequency, the discretization in the 
TDM is not compensated for in the HBM calculations. 

Then there should be agreement between the sampling 
frequency f s in the TDM and the number of samples N t 
per period in the HBM. Their relation is given by 



N t = 



(33) 



In order to have a sufficiently high sampling rate, we have 
chosen N t = 512. The playing frequency f p is plotted in 
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FIG. 14. Comparison between HBM and TDM of the am- 
plitude of (a) the first harmonic Pi and (b) the dimensionless 
playing frequency f P /f r as the blowing pressure 7 increases for 
a clarinet-like system with viscous losses and weak and normal 
reed interaction. TDM values for £ == 0.50 and 7 > 0.48 are 
omitted due to period doubling. So are the beating regimes of 
HBM calculations. (/„ = 51100Hz, 7Vt = 512, f r = 103.4 Hz, 
N p = 15) 
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FIG. 15. The pressure (a) and volume-flow (b) wave form of 
the Helmholtz solution and a 3-level sister solution calculated 
by the HBM employing the simple clarinet model with £ = 
0.5, 7 = 0.4, N p = 99, 77 = 10 -5 . 



Figure fHb. and we used an average f s = 51100 Hz for 
both the HBM and the TDM. 

Finally, it seems also necessary that N p and N t are 
chosen so that 

N p + l = -£. (34) 

In practice, however, when comparing bifurcation dia- 
grams of the first harmonic Pi , as in Figure El rather 
low values of N p give good results. Nevertheless, more 
harmonics are obviously needed to compare waveforms 
in the time domain, especially far from the oscillation 
threshold. 

V. PRACTICAL EXPERIENCES 

A. Multiple solutions 

As we consider a nonlinear problem, we cannot antici- 
pate the number of solutions. Therefore, it should not be 



surprising that it is possible to obtain multiple solutions 
for a given set of parameter values. When searching for 
a particular solution, this may be a practical problem. 
Fritz et al. ? have discovered that some solutions seem 
to disappear when increasing the number of harmonics 
N p , implying that solutions may arise from the trunca- 
tion to a finite N p . We have now discovered alternative 
solutions that persist even at very high N p . 

Let us illustrate this with the simple model of the clar- 
inet used in Section II V H II where the reed is a spring 
without mass or damping, the nonlinearity is given by 
equation (|23fl . and the bore is an ideal cylinder with 
nearly lossless propagation and no dispersion. Figure H5l 
shows a three-level sister solution together with the re- 
lated Helmholtz solution for a large number of harmonics, 
N p = 2000. 

A solution of the lossless problem should satisfy the 
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criteria - 



r P (t + tt) = - P (t) 

\ U(t + 7r) = lt(t) 



(35) 



(the dimcnsionless period being 2tt), as well as the 
conditions stated before equation (|24J) . noting that 
p(t) = u{t) = for all t is the static solution. It is eas- 
ily verified graphically that both of the solutions in Fig- 
urc^]satisfy these conditions. Moreover, since they also 
satisfy equation 1)23(1 , the three- level solution is a solution 
of the lossless model. 

Whereas the system of time-domain equations (|35|) has 
an infinity of solutions, truncation in frequency-domain 
limits the number of solutions. The unique solution of 
the HBM with only one harmonic is obviously a sine. 
Let us analyse the situation in the simplest nontrivial 
case of the lossless problem with two odd harmonics and 
a cubic expansion for nonlinear coupling. Ignoring even 
harmonics, the HBM gives a system of two equations (see 
Kergomard et al.) ? : 



a = 3Pi2(l 
ax = Pi2(H 



f.x + 2|x| 2 ) (a) 
3x|x| 2 + 6x), (b) 



(36) 



where a = —A/C and x = Ps/Pi. As equation Q36fa ) 
imposes P3 to be real, solving this system amounts to 
solving 



x 3 + x 2 



1/3. 



(37) 



This equation has three real solutions x ~ —1.5151, 
-0.2776 and 0.7926. All of them are found by Harm- 
bal for negligible losses (77 = 10~ 5 ), and the correspond- 
ing waveforms are presented in Figure El We note that 
the second solution leads to the Helmholtz motion when 
increasing the number of harmonics (with the theoret- 
ical value known to be x = —1/3) whereas the third 
one corresponds to the three-level solution in Figure El 
We can also easily imagine that these three solutions of 
the truncated problem are three-harmonic approxima- 
tions of square waves that are distributed on three levels: 
~ ±0.5 and p = 0. Respectively, they have two, zero, 
and one steps at the zero-level. It should be noted that 
the conditions for the continuous problem do not 
constrain the duration of each step. Figure IT71 shows two 
such twin solutions for N p = 99 corresponding to the 
three-level solution in Figure El This has to be kept in 
mind when increasing N p using the HBM. 

While the Helmholtz motion is known to be stable, 7 
the two three-level solutions can be considered as a com- 
bination of the static solution (the zero level) and the 
square wave (two levels with opposite values). Since we 
know from Kergomard 7 that in the case of ideal propaga- 
tion (neither losses nor dispersion), the stability domain 
of these two solutions are mutually exclusive, it can be 
concluded that the three-level solutions are unstable. 

Taking into account losses in the propagation does not 
make the three-level solutions vanish. But a simple rea- 
soning to determine the stability of this solution is not 
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FIG. 16. The pressure waveform of the three solutions found 
by the HBM with N p = 3 employing the simple clarinet model 
with C = 0.5, 7 = 0.4, r) = 10" 5 . 
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FIG. 17. The pressure waveform of two solutions that differ 
by the duration of their steps, found by the HBM employing 
the simple clarinet model with £ = 0.5, 7 = 0.4, r\ = 10 _J , 
N p = 99. 



possible in this case. To the authors knowledge, however, 
such a solution has never been observed experimentally 
at low level of excitation. 



B. Initial value of the playing frequency 

A practical difficulty encountered is the convergence of 
the playing frequency f p . If its initial value is not close 
enough to the solution, divergence is almost inevitable. 
This occurs because the resonator impedance Z r tends 
to vanish outside the immediate surroundings of the res- 
onance peaks of the resonator, rendering F(P, f p ) very 
small and thereby G ~ PjP\ nearly constant with re- 
spect to f p . The slope dG/df p thus becomes close to 
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zero, the Newton step leads far away from the solution, 
and convergence fails. Dissipation widens the resonance 
peaks and thus also the convergence range. 

For a simple system where the playing frequency is 
known to correspond to a resonance peak of the tube, 
initializing f p is easy. However, with dispersion or other 
inharmonic effects, choosing an initial value for f p may 
be difficult. In Harmbal the problem may to some extent 
be avoided by the possibility of gradually adding the dis- 
persion (or other inharmonic effects), so that the play- 
ing frequency can be followed quasi-continuously from a 
known solution without dispersion, for instance by using 
hbmap. 

VI. CONCLUSIONS 

The harmonic balance method (HBM) is suited for 
studies of self-sustained oscillations of musical instru- 
ments, and the computer program Harmbal has been de- 
veloped for this application. It is available with its source 
code, 7 has a free licence, and is already in use by sev- 
eral researchers. It is programmed in C, runs fast, and is 
easily used by other application, such as for continuation 
purposes. 

Some difficulties are related to the digital sampling of 
the signal and can be solved by introducing a backtrack- 
ing mechanism. When using a large number of harmon- 
ics, the extreme case of the (lossless) Helmholtz motion 
can be solved for different shapes of resonators. Never- 
theless, the value of the first harmonic Pi seems to be 
well predicted by lower values of N p , in particular close 
to the threshold of a direct bifurcation. For the saxo- 
phone we used a stepped-cone bore and observed one or 
more dips during the longest part of the period, depend- 



ing on the number of steps. These dips approach pure 
impulses as N p increases. The number of samples N t in 
a period showed to be insignificant for these dips. 

The HBM can lead to some alternative solutions for a 
unique set of parameters. The nondissipative versions of 
these solutions satisfy the continuous model equations, 
but they are not stable and thus cannot be attained by 
ab initio time-domain calculations. Another problem is 
the great sensitivity to the guessed playing frequency. 

As a consequence, a certain expertise is needed in order 
to use the method, but, thanks to an automatic contin- 
uation procedure, the calculation is easy. We note that 
also experimental results can be used for the impedance 
of the resonator. 
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1 Self-sustained is a term indicating oscillation driven by a constant 
energy input. 

2 The term non-interactive means that the user has no influence on 
the program while it is running. 
3 This result corrects equation (14) in ref.? 
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